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ABSTRACT 

Two-dimensional  simulations  of  the  1 1  January  1972  Boulder,  Colorado,  windstorm,  obtained  from  1 1  diverse 
nonhydrostatic  models,  are  intercompared  with  special  emphasis  on  the  turbulent  breakdown  of  topographically 
forced  gravity  waves,  as  part  of  the  preparation  for  the  Mesoscale  Alpine  Programme  field  phase.  The  sounding 
used  to  initialize  the  models  is  more  representative  of  the  actual  lower  stratosphere  than  those  applied  in  previous 
simulations.  Upper-level  breaking  is  predicted  by  all  models  in  comparable  horizontal  locations  and  vertical 
layers,  which  suggests  that  gravity  wave  breaking  may  be  quite  predictable  in  some  circumstances.  Characteristics 
of  the  breaking  include  the  following:  pronounced  turbulence  in  the  13-16-km  and  18-20-km  layers  positioned 
beneath  a  critical  level  near  21 -km,  a  well-defined  upstream  tilt  with  height,  and  enhancement  of  upper-level 
breaking  superpositioned  above  the  low-level  hydraulic  jump.  Sensitivity  experiments  indicate  that  the  structure 
of  the  wave  breaking  was  impacted  by  the  numerical  dissipation,  numerical  representation  of  the  horizontal 
advection,  and  lateral  boundary  conditions.  Small  vertical  wavelength  variations  in  the  shear  and  stability  above 
10  km  contributed  to  significant  changes  in  the  structures  associated  with  wave  breaking.  Simulation  of  this 
case  is  ideal  for  testing  and  evaluation  of  mesoscale  numerical  models  and  numerical  algorithms  because  of  the 
complex  wave-breaking  response. 


1.  Introduction 

The  basic  theory  for  inviscid  mountain  waves  forced 
by  airflow  over  a  two-dimensional  obstacle  in  a  stratified 
atmosphere  has  been  established  for  several  decades  and 
has  been  the  subject  of  several  reviews  (e.g.,  Queney 
et  al.  1960;  Smith  1979).  As  vertically  propagating  in¬ 
ternal  waves  amplify,  in  part  due  to  the  decrease  of 
atmospheric  density  with  altitude  and  nonlinear  pro¬ 
cesses,  the  mountain  waves  may  overturn  and  lose  en¬ 
ergy  due  to  turbulent  breakdown.  Although  two-  and 
three-dimensional  models  have  been  able  to  successfully 
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simulate  gravity-wave  generation  (e.g..  Smith  1979), 
evolution  (e.g..  Smith  and  Grpnas  1993;  Schar  and  Out¬ 
ran  1997),  and  breakdown  (Clark  and  Peltier  1977;  Bac- 
meister  and  Schoeberl  1989),  numerous  questions  exist 
regarding  the  turbulent  breakdown  of  three-dimensional 
gravity  waves  in  nature.  Furthermore,  it  has  not  been 
conclusively  demonstrated  that  numerical  models  can 
accurately  predict  the  occurrence  of  wave  breaking  over 
three-dimensional  complex  topography. 

Mountain-wave  breaking  has  an  important  influence 
on  the  atmosphere  for  several  reasons.  The  vertical  mo¬ 
mentum  flux  associated  with  wave  breaking  and  oro¬ 
graphic  drag  has  a  profound  impact  on  the  atmospheric 
general  circulation  (Bretherton  1969)  and  requires  care¬ 
ful  representation  in  general  circulation  and  numerical 
weather  prediction  models  (Palmer  et  al.  1986;  Lott 
1995).  Downslope  windstorms  and  drag  enhancement 
may  occur  as  a  result  of  resonant  mountain-wave  am- 
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plification  positioned  beneath  upper-level  wave-break¬ 
ing  regions  that  contain  “wave-induced”  critical  levels 
(Clark  and  Peltier  1977;  Peltier  and  Clark  1979).  Peltier 
and  Clark  (1979)  attribute  the  windstorm  formation  to 
the  bifurcation  character  of  the  flow  response  governed 
by  the  occurrence  of  breaking.  Large  vertical  diffusivity 
associated  with  the  turbulent  breakdown  of  orographi- 
cally  generated  gravity  waves  results  in  efficient  vertical 
mixing  of  water  vapor,  aerosols,  various  chemical  spe¬ 
cies,  and  trace  gasses  (Dornbrack  and  Diirbeck  1998). 
Wave  breaking  represents  an  internal  dissipation  mech¬ 
anism  that  leads  to  the  generation  of  potential  vorticity 
(Schar  and  Smith  1993)  and  subsequently  may  influence 
the  downstream  mesoscale  and  synoptic-scale  flow 
(Thorpe  et  al.  1993;  Aebischer  and  Schar  1998).  Clear- 
air  turbulence  in  conjunction  with  the  breakdown  of 
gravity  waves  occurs  frequently  over  complex  terrain 
(Nastrom  and  Fritts  1992)  and  is  a  significant  aviation 
hazard  (Lilly  1978;  Ralph  et  al.  1997). 

One  of  the  goals  of  the  Mesoscale  Alpine  Programme 
(MAP)  is  to  obtain  a  more  complete  understanding  of 
the  dynamics  governing  the  turbulent  breakdown  of  or- 
ographically  generated  gravity  waves.  MAP  is  a  coor¬ 
dinated  international  effort  to  examine  moist  and  dry 
processes  associated  with  three-dimensional  complex- 
topographic  flows  including  the  physics  and  dynamics 
of  heavy  precipitation  events,  potential  vorticity  gen¬ 
eration,  gap  flow,  as  well  as  gravity  wave  breaking 
(Binder  and  Schar  1996;  Houze  et  al.  1998).  High-res- 
olution  numerical  simulations  (Ax  ~  3-6  km  and  ~50 
vertical  levels)  were  conducted  on  an  operational  basis 
during  the  MAP  intensive  observation  periods  in  an 
attempt  to  predict  the  altitude  and  location  of  wave 
breaking  regions  over  the  Alps.  The  numerical  forecasts 
were  used  as  guidance  for  research  aircraft  missions 
designed  to  observe  mountain  wave  breaking. 

In  preparation  for  the  MAP  held  phase,  a  collabo¬ 
rative  effort  was  established,  as  part  of  the  pre-MAP 
research  effort,  to  assess  the  ability  of  high-resolution 
numerical  models  to  predict  the  location  and  strength 
of  breaking  waves  in  a  series  of  experiments  involving 
airflow  over  topography.  The  first  test  is  a  multimodel 
intercomparison  of  two-dimensional,  high-resolution 
numerical  simulations  of  the  11  January  1972  down- 
slope  windstorm  that  occurred  along  the  Colorado  Front 
Range.  The  dynamics  of  this  windstorm  has  been  the 
subject  of  numerous  theoretical  and  numerical  studies 
(e.g.,  Klemp  and  Lilly  1975;  Clark  and  Peltier  1977; 
Klemp  and  Lilly  1978;  Peltier  and  Clark  1979;  Peltier 
and  Clark  1983;  Clark  and  Farley  1984;  Durran  1986). 
The  primary  focus  of  this  study  is  on  the  predictability 
and  characteristics  of  the  numerically  simulated  wave 
breaking  rather  than  the  dynamics  of  the  leeside  wind¬ 
storm.  Multiple  regions  of  complex  upper-level  wave 
breaking  develop  in  this  case  and  represent  a  challeng¬ 
ing  test  of  the  capability  of  various  model  architectures 
to  produce  equivalent  wave-breaking  simulations  using 
identical  initial  conditions.  The  goal  of  this  work  is  to 


examine  the  variability  in  the  simulated  gravity  wave 
breaking  among  the  1 1  numerical  models  applied  to  this 
case  rather  than  focusing  on  replication  the  flnescale 
wind  storm  aspects.  The  distribution  of  the  wave-break¬ 
ing  regions  is  particularly  complex  in  the  stratosphere, 
and  it  follows  that  one  focus  will  be  on  the  lower-strato¬ 
spheric  flow.  In  order  to  faithfully  simulate  the  strato¬ 
spheric  response,  the  models  were  initialized  using  ob¬ 
servations  from  the  upstream  sounding  (Grand  Junction, 
Colorado)  through  an  altitude  of  25  km.  As  will  be 
discussed,  this  initial  sounding  is  much  more  represen¬ 
tative  of  the  actual  lower  stratospheric  structure  than 
those  used  in  previous  simulations  of  the  11  January 
1972  windstorm. 

2.  Experimental  design  and  description  of 

numerical  models 

The  configuration  of  the  models  is  standardized  as 
much  as  possible  to  facilitate  a  meaningful  intercom¬ 
parison.  The  models  are  applied  in  a  two-dimensional 
mode  with  a  horizontal  grid  increment  of  1  km  and  a 
vertical  grid  spacing  of  200  m  up  to  an  altitude  of  25 
km  (~125  vertical  levels).  The  mountain  profile  is  a 
witch  of  Agnesi  curve. 


where  h  is  the  mountain  height  and  a  is  the  mountain 
half  width.  In  these  simulations,  a  =  10  km  and  h  = 
2  km,  which  are  typical  values  for  the  Colorado  Front 
Range.  Although  portions  of  the  continental  divide  up¬ 
stream  of  Boulder  can  be  considered  quasi-two-dimen- 
sional,  the  numerical  models  are  applied  in  a  two-di¬ 
mensional  mode  in  this  study  to  facilitate  simple  and 
fundamental  comparisons  of  the  wave-breaking  simu¬ 
lations.  The  horizontal  boundaries  are  located  10a  (100 
km)  upstream  and  12a  (120  km)  downstream  of  the 
mountain.  Moist  effects  and  rotation  are  not  considered. 
A  free-slip  lower  boundary  condition  is  used.  The  model 
results  are  compared  at  the  3-h  time  of  a4-h  integration. 
The  results  at  the  3-h  time  are  characterized  by  varia¬ 
tions  on  slow  timescales,  with  all  models  obtaining  nom¬ 
inally  equivalent  predictions  of  the  windstorm  and  up¬ 
per-level  wave  breaking. 

The  initial  conditions  are  horizontally  homogeneous 
and  based  upon  the  upstream  1200  UTC  11  January 
1972  Grand  Junction,  Colorado,  sounding  shown  in 
Figs.  la,b.  The  sounding  contains  multiple  shear  layers 
(Fig.  la),  strong  low-level  cross-mountain  winds  (~15 
m  s  ')>  and  a  well-defined  stable  layer  above  the  moun¬ 
tain  crest  between  2  km  and  7  km  (Fig.  lb),  with  the 
latter  two  characteristics  favorable  for  downslope  wind¬ 
storms  (e.g.,  Brinkman  1974;  Durran  1986).  A  critical 
level  is  present  at  ~21  km  where  the  cross-mountain 
wind  speed  is  small.  It  is  noteworthy  that  the  wind  and 
stability  structure  of  the  original  sounding  is  retained 
up  to  25  km,  in  contrast  to  previous  numerical  studies 
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Fig.  1.  Vertical  profiles  of  the  (a)  w-component  (m  s“')  and  (b)  temperature  (K)  from  the  Grand 
Junction,  CO,  sounding  for  1200  UTC  11  January  1972  (heavy  solid  line)  and  initial  condition 
sensitivity  experiment  (dotted  lines).  Vertical  profiles  of  the  (c)  w-component  (m  s“')  and  (d) 
temperature  (K)  used  in  previous  simulations  of  the  1200  UTC  11  January  1972  windstorm  (after 
Peltier  and  Clark  1979). 


that  typically  used  smoothed  profiles  in  the  stratosphere 
(e.g.,  Peltier  and  Clark  1979),  as  shown  in  Figs.  lc,d. 
The  initial  sounding  applied  in  the  present  study  con¬ 
tains  a  more  representative  initial  state  in  the  strato¬ 
sphere  than  used  in  previous  simulations  of  the  1 1  Jan¬ 
uary  1972  windstorm,  which  should  presumably  lead  to 
a  more  realistic  simulation  of  upper-level  wave  break¬ 
ing. 

The  11  January  1972  windstorm  along  the  Colorado 
Front  Range  has  been  analyzed  on  numerous  occasions 
both  observationally  (Lilly  and  Zipser  1972;  Lilly  1978) 
and  numerically  (beginning  with  Klemp  and  Lilly  1975; 
Klemp  and  Lilly  1978;  Peltier  and  Clark  1979)  because 
of  the  rare  direct  in  situ  measurements  of  the  mountain 
wave  structure.  The  composite  analyses  of  potential 
temperature  and  westerly  wind  component  obtained 
from  Sabreliner  and  Queen- Air  aircraft  flight  data  for 
this  situation  are  shown  in  Fig.  2  (Klemp  and  Lilly 


1975).  A  large-amplitude  wave  extends  over  300  hPa 
in  the  upper  troposphere  (Fig.  2a)  and  is  nearly  coin¬ 
cident  with  a  deep  layer  of  weak  or  near-reversed  cross¬ 
mountain  flow  (Fig.  2b),  which  is  likely  a  manifestation 
of  mountain-wave  amplification  and  breaking.  Lilly 
(1978)  surmises  that  the  60  m  s  '  downslope  winds  in 
the  lee  are  consistent  with  a  deep  layer  between  200 
and  700  hPa  upstream  of  the  mountains  passing  through 
a  260-hPa  layer  immediately  above  the  surface  in  the 
lee  of  the  mountain  crest. 

The  numerical  model  characteristics  and  architectures 
applied  in  this  study  vary  considerably,  as  summarized 
in  Table  1.  The  models  used  in  the  intercomparison 
include  the  following;  the  Advanced  Regional  Predic¬ 
tion  System  (ARPS)  (Xue  et  al.  1995),  Coupled  Ocean/ 
Atmosphere  Mesoscale  Prediction  System  (COAMPS) 
(Hodur  1997),  Durran  and  Klemp  (1983)  model  (DK83), 
Eulerian/semi-Lagrangian  model  (EULAG)  (Smolarkie- 
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Fig.  2.  Analysis  of  (a)  potential  temperature  and  (b)  w-component  speed  from  aircraft  flight  data  and  sondes  taken  on  11  January  1972. 
Aircraft  tracks  are  shown  by  dashed  lines  with  locations  of  significant  turbulence  shown  by  the  plus  signs.  The  heavy  dashed  line  separates 
data  taken  by  the  Queen  Air  (before  2200  UTC)  from  the  Saberliner  data  (after  0000  UTC)  (from  Klemp  and  Lilly  1975). 
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Table  1.  Model  formulation  and  numerical  characteristics.  The  model  equations  sets  are  delineated  either  as  compressible  (C),  anelastic 
(A),  or  Arakawa  Jacobian  (A-J).  The  finite  difference  algorithms  are  referred  to  as  centered  in  time  (CIT),  forward  in  time  (FIT),  centered 
in  space  (CIS),  upstream  in  space  (UIS),  and  semi-Lagrangian  (SL). 


Model 

Eqn. 

set 

Vert. 

coord. 

Time 

differ. 

Time  step/ 
filter 
coef.  (y) 

Horizon. 

advect. 

Horizon. 

smooth/ 

Ko 

Top 

bound. 

cond. 

Initial. 

Turb. 

param. 

ARPS 

c 

CIT 

Large:  9  s 

4th  Order 

4th  Order 

Sponge 

1.5  Order 

Semi-implicit 

Small:  1  s 

CIS 

(TKE) 

y  =  0.05 

1.1  X  10*  m*  s-' 

COAMPS 

c 

CIT 

Large:  3  s 

2d  Order 

4th  Order 

Sponge 

g,  U 

1.5  Order 

Semi-implicit 

Small:  1.5  s 

CIS 

Spin  Up 

(TKE) 

y  =  0.2 

1.3  X  10'>  m-*  s-' 

CUMM 

A 

(7. 

LIT 

Large:  1  s 

7  th  Order 

Implicit/ 

Sponge 

1.5  Order 

Small:  N/A 

SL 

Shapiro 

(TKE) 

N/A 

Eilter 

DK83 

c 

cr. 

CIT 

Large:  12  s 

4th  Order 

4th  Order 

Radiative 

g,  u 

1st  Order 

Simi-implicit 

Small:  3  s 

CIS 

Spin  Up 

y  =  0.2 

1.3  X  10’ m*  s-‘ 

EULAG 

A 

LIT 

Large:  5  s 

2d  Order 

Implicit/2d 

Sponge 

Potential 

1.5  Order 

Small:  N/A 

UIS 

Order 

flow 

(TKE) 

N/A 

8  m"  s“' 

(2  D) 

MESO-NH 

A 

CIT 

Large:  8  s 

2d  Order 

4th  Order 

Sponge 

1.5  Order 

Small:  N/A 

CIS 

(TKE) 

y  =  0.2 

4.1  X  10*  m*  s-' 

MM5 

C 

O', 

CIT 

Large:  3  s 

2d  Order 

4th  Order 

Radiative 

1st  Order 

Semi-implicit 

Small:  1  s 

CIS 

y  =  0.1 

1.0  X  10’ m’ s"' 

NTU/Purdue 

c 

2  Step 

Large:  1  s 

2d  Order 

2d  Order 

Sponge 

g,  u 

1st  Order 

Forward-Backward 

Small:  0.3  s 

CIS 

Spin  Up 

N/A 

1.3  X  10’ m’ s-' 

RAMS 

c 

O', 

CIT 

Large:  3  s 

2d  Order 

4th  Order 

Sponge 

U 

1.5  Order 

Semi-implicit 

Small:  1  s 

CIS 

Spin  Up 

(TKE) 

y  =  0.2 

RIMS 

c 

Adams-Bash: 

Large:  1  s 

4th  Order 

4th  Order 

Sponge 

1st  Order 

Semi-implicit 

Small:  0.2  s 

CIS 

y  =  N/A 

5.0  X  10’ m’ s-' 

UCLA 

A-J 

step 

CIT 

Large:  3  s 

2d  Order 

None 

Sponge 

None 

Small:  N/A 

CIS 

y  =  N/A 

wicz  and  Margolin  1997),  the  Central  University  Me- 
soscale  Model  (CUMM)  (Huang  2000),  Mesoscale- 
Nonhydrostatic  (MESO-NH)  Model  (Lafore  et  al. 
1998),  Penn  State-NCAR  Mesoscale  Model  (MM5) 
(Dudhia  1993),  National  Taiwan  University  (NTU)/Pur- 
due  Model  (Hsu  and  Sun  1999,  manuscript  submitted 
to  Mon.  Wea.  Rev.),  Regional  Atmospheric  Modeling 
System  (RAMS)  (Pielke  et  al.  1992),  Regional  Inte¬ 
grated  Modeling  System  (RIMS)  (Chen  1991),  and  the 
UCLA  model  (Landau  et  al.  1997).  All  of  the  models 
applied  in  this  study  are  nonhydrostatic  with  several 
making  use  of  the  anelastic  equation  set  (EULAG, 
CUMM,  MESO-NH),  others  solving  the  fully  com¬ 
pressible  equations  (ARPS,  COAMPS,  DK83,  MM5, 
NTU/Purdue,  RAMS,  RIMS)  and  one  model  applying 
the  Arakawa  Jacobian  method  to  predict  the  stream- 
function  and  vorticity  (UCLA).  The  vertical  coordinates 
are  terrain  following  (cr.)  with  the  exception  of  the 
UCLA  model,  which  uses  a  block-type  coordinate  sys¬ 
tem.  The  models  differ  significantly  in  spatial  and  tem¬ 
poral  finite  differencing  techniques.  Lor  example,  a 
group  of  models  (COAMPS,  MM5,  RAMS)  make  use 
of  standard  formulations  such  as  second-order  accurate 


centered  temporal  and  spatial  differencing,  in  contrast 
to  other  models  that  use  forward-temporal  and  up- 
stream-spatial  differencing  (EULAG)  or  centered  tem¬ 
poral  and  higher-order  accurate  spatial  differencing 
(ARPS,  CUMM,  DK83).  Significant  differences  among 
models  exist  in  the  magnitude  of  the  horizontal  smooth¬ 
ing  and  time  filter  coefficients.  Additionally,  some  mod¬ 
els  employ  second-order  horizontal  smoothing  (EU¬ 
LAG,  NTU/Purdue)  while  others  use  fourth-order 
smoothing.  Turbulence  is  represented  by  simple  first- 
order  closure  techniques  in  several  models  (DK83, 
MM5,  RIMS,  NTU/Purdue).  In  contrast,  an  explicit  pre¬ 
dictive  equation  for  turbulent  kinetic  energy  is  employed 
by  a  number  of  models  (ARPS,  COAMPS,  CUMM, 
EULAG,  MESO-NH,  RAMS)  to  parameterize  the  sub- 
gridscale  vertical  transport.  It  is  noteworthy  that  the 
UCLA  model  does  not  parameterize  the  subgridscale 
turbulence  and  relies  on  the  gridscale  diffusive  pro¬ 
cesses  for  vertical  diffusion  and  mixing.  Several  models 
use  a  spinup  period  for  the  M-component  and/or  gravity 
(CO AMPS,  DK83,  NTU/Purdue,  RAMS),  while  EU¬ 
LAG  uses  a  two-dimensional  potential  flow  initializa¬ 
tion.  The  UCLA  model  makes  use  of  cyclic  lateral 
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boundary  conditions  in  contrast  to  the  other  models, 
which  use  radiative  boundary  conditions. 

3.  Results 

a.  Control  simulations 

The  simulated  potential  temperature  after  3  h  for  the 
1 1  models  is  shown  in  Fig.  3.  All  of  the  models  produce 
two  deep  regions  of  wave  breaking  located  in  the  strato¬ 
sphere,  near  —13-16  km  and  —18-20  km,  with  a  ver¬ 
tical  wavelength  of  —3-4  km  and  an  upstream  vertical 
tilt.  The  UCLA  model  results  (Fig.  3k)  differ  from  the 
other  simulations  in  that  the  stratospheric  wave-break¬ 
ing  regions  are  less  coherent,  with  numerous  small-scale 
superadiabatic  regions,  particularly  in  the  troposphere 
(some  hltering  has  been  applied  to  the  UCLA  model 
output).  The  lack  of  a  subgridscale  eddy  viscosity  pa¬ 
rameterization  may  lead  to  greater  high-frequency  struc¬ 
tures  in  the  UCLA  model  relative  to  the  other  simula¬ 
tions,  which  make  use  of  subgridscale  turbulence 
schemes  (Table  1). 

All  1 1  of  the  models  produce  a  signihcant  leeside 
windstorm,  as  indicated  in  Fig.  4.  The  lee-slope  wind 
maximum  at  the  3-h  time  varies  considerably  from  —86 
m  s  '  for  RIMS  to  —43  m  s  '  for  the  UCLA  model. 
This  variability  in  the  model  simulations  may  be  a 
reflection  of  differences  in  model  formulation  as  well 
as  the  unsteady  character  of  the  wave  breaking.  Sig¬ 
nificant  differences  among  the  model  simulations  exist 
in  the  flow  upstream  of  the  barrier  near  the  7  km  al¬ 
titude,  ranging  from  —43  m  s^'  for  the  MM5  results 
to  greater  than  52  m  s  '  for  CUMM.  A  well-defined 
hydraulic  jump,  as  evidenced  in  the  potential  temper¬ 
ature  and  velocity  fields,  is  positioned  along  the  lee 
slope  in  several  model  simulations  (e.g.,  DK83,  EU- 
LAG)  and  10-25  km  downstream  in  other  simulations 
(CUMM,  RAMS,  RIMS).  Substantial  variations  exist 
in  the  movement  of  the  jump  and  the  transition  to  the 
ambient  downstream  conditions.  These  differences  are 
consistent  with  the  large  model  variability  (factor  of 
2)  in  maximum  vertical  velocity  near  the  hydraulic 
jump  at  the  3-h  time.  Variations  in  the  hydraulic  jump 
behavior  and  characteristics  are  likely  influenced  by 
lateral  boundary  formulation  and  the  unsteady  nature 
of  the  solution. 

The  model  simulations  reproduced  the  characteristics 
of  the  tropospheric  structure  analyzed  by  Klemp  and 
Lilly  (1975)  and  Lilly  (1978)  (Fig.  2)  reasonably  well, 
such  as  the  wind  speed  maximum  in  the  lee,  the  hy¬ 
draulic  jump,  the  midtropospheric  cross-mountain  wind 
component  minimum  in  the  lee,  and  the  weak  stability 
layer  associated  with  the  midlevel  gravity  wave  steep¬ 
ening  and  breaking.  However,  Klemp  and  Lilly’s  anal¬ 
ysis  does  not  support  the  model  simulated  downward 
tropopause  depression  associated  with  breaking  super- 
positioned  over  the  low-level  hydraulic  jump.  The  hy¬ 
draulic  jump  extends  to  a  greater  altitude  and  has  a 


steeper  slope  in  nearly  all  of  the  simulations  than  in¬ 
dicated  by  the  aircraft  observations.  Some  of  these  dif¬ 
ferences  may  be  related  to  the  simplicity  of  model  ar¬ 
chitectures,  including  the  two-dimensional  and  inviscid 
configuration,  as  well  as  the  simplified  model  initial 
state  and  topography  representation. 

The  upper-level  wave  breaking  is  apparent  in  the 
Richardson  number  (R-)  fields  in  Fig.  5  [sgn(R,)IR,l''  ^ 
was  chosen  for  graphical  purposes].  The  upper-level 
breaking  regions  in  several  of  the  model  simulations 
extend  to  —  lOa  downstream  and  —2a  upstream  of  the 
mountain.  In  all  of  the  simulations,  with  the  exception 
of  the  UCLA  model,  the  largest  amplitude  breaking,  as 
measured  by  the  vertical  excursion  of  the  isentropes 
(Fig.  3),  is  located  just  downstream  and  above  the  lower- 
tropospheric  hydraulic  jump.  Regions  of  small  R,  can 
be  found  in  the  lower  troposphere  in  the  vicinity  of  the 
hydraulic  jump  and  above  the  lee  slopes  in  a  nearly 
adiabatic  layer  between  6  and  8  km  (Fig.  3).  This  layer 
is  characterized  by  flow  stagnation,  with  some  models 
containing  a  deep  region  of  reversed  flow  (Fig.  4).  The 
damping  effect  of  the  critical  layer,  which  limits  the 
vertical  extent  of  the  breaking,  is  apparent  in  the  R,  (Fig. 
5)  and  potential  temperature  fields  (Fig.  3).  In  spite  of 
significant  differences  between  the  models  in  the  low- 
level  hydraulic  structures,  vertical  velocities,  and  leeside 
winds,  the  vertical  and  horizontal  locations  of  the  wave 
breaking  are  generally  similar  in  all  model  simulations. 
The  UCLA  model,  which  produced  breaking  structures 
in  similar  layers  as  the  other  models,  however,  differed 
substantially  in  terms  of  the  vertical  coherence  and  high- 
frequency  waves.  Some  of  the  differences  between  mod¬ 
el  simulations  may  arise  due  to  the  transient  nature  of 
the  wave  breaking,  which  makes  comparison  at  a  par¬ 
ticular  time  difficult.  Even  though  a  complete  steady 
state  is  generally  not  achievable  in  these  simulations, 
comparison  of  the  model  solutions  at  a  given  time  is 
still  meaningful,  especially  with  regard  to  predictability 
issues. 

Although  there  are  many  broad  similarities  in  the 
overall  structure  of  the  wave  breaking  in  nearly  all  of 
the  simulations,  some  noteworthy  differences  in  the  up¬ 
per-level  breaking  structures  exist.  Eor  example,  several 
model  simulations  contain  breaking  regions  in  the  lower 
stratosphere  downstream  of  the  mountain  that  penetrate 
downward  below  7  km  resulting  in  a  deep  adiabatic 
layer  (Eigs.  3a  and  3j).  Additionally,  these  two  simu¬ 
lations  contain  more  small-scale  structures  associated 
with  the  wave  breaking,  which  may  be  due  to  the 
strength  of  the  numerical  smoothing.  A  number  of  mod¬ 
el  simulations  contain  trapped  waves  downstream  of  the 
wave  breaking  in  the  10-14  km  layer  (Eigs.  3e,  3g,  3h, 
and  3i).  Some  of  the  differences  between  model  solu¬ 
tions  may  be  related  to  the  relative  strength  of  the  ex¬ 
plicit  or  implicit  numerical  smoothing  (Table  1).  Some 
models  that  used  smaller  horizontal  smoothing  eddy  dif- 
fusivity  (Kfl)  and  Robert  (1966)/Asselin  (1972)  time 
Alter  (y)  coefficients  tend  to  contain  smaller-scale  wave 
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Eig.  3.  Vertical  cross  section  of  the  simulated  potential  temperature  after  3  h  for  (a)  ARPS,  (b)  COAMPS,  (c)  CUMM,  (d)  DK83,  (e)  EULAG,  (f)  MESO-NH,  (g)  MM5,  (h)  NTU/Purdue, 

(i)  RAMS,  (j)  RIMS,  and  (k)  UCLA  models.  The  contour  interval  is  8  K. 
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RAMS,  (j)  RIMS,  and  (k)  UCLA  models.  The  contour  interval  is  8  m  s 
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Eig.  5.  The  sgn(R,.)IR,.|“=  for  the  3-h  time  for  (a)  COAMPS,  (b)  DK83  model,  (c)  EULAG,  and  (d)  RAMS.  The  isopleths 

between  1.0  and  —5  are  displayed  every  0.5. 


breaking  substructures  relative  to  the  other  models  that 
apply  a  larger  and/or  y  (see  Figs.  3,  4).  However, 
there  are  notable  exceptions  to  this  relationship  clearly 
suggesting  aspects  of  the  model  formulation  other  than 
horizontal  smoothing  likely  control  the  simulation  of 
hnescale  breaking  signatures.  For  example,  RIMS  (Figs. 
3j  and  4j),  which  used  a  relatively  large  qualitatively 
produces  higher  amplitude  small-wavelength  structures 
than  MESO-NH  (Figs.  3f  and  4f),  which  makes  use  of 
a  relative  small  K^. 

b.  Sensitivity  tests 

Experiments  were  performed  to  isolate  the  sensitivity 
of  the  wave  breaking  to  the  dissipation,  advection 
scheme,  and  lateral  boundary  conditions.  The  3-h  po¬ 
tential  temperature  held,  obtained  using  CO  AMPS  with 


second-order  horizontal  smoothing,  is  shown  in  Eig.  6a. 
In  comparison  with  the  control  simulation  (Eig.  3b), 
which  made  use  of  a  more  scale-selective  fourth-order 
horizontal  smoothing,  not  surprisingly  the  increased  dis¬ 
sipation  leads  to  shallower  neutral  layers  downstream 
of  the  major  wave-breaking  source  regions.  The  vertical 
excursion  of  the  isentropes  downstream  of  the  hydraulic 
jump  is  reduced  as  well.  Moreover,  the  sensitivity  of 
these  breaking  structures  to  the  horizontal  smoothing  or 
diffusion  is  illustrated  by  the  ARPS  results  in  Pig.  6b, 
which  is  a  simulation  that  made  use  of  a  fourth-order 
horizontal  smoothing  coefficient  that  was  doubled  rel¬ 
ative  to  the  control  simulation  (Fig.  3a).  As  expected, 
the  increase  in  the  horizontal  smoothing  results  in  a 
decrease  of  the  small-scale  structure  associated  with  the 
breaking  and  a  decrease  in  the  amplitude  of  the  breaking 
positioned  above  the  hydraulic  jump.  An  increase  in  the 
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Fig.  6.  The  3-h  potential  temperature  for  (a)  COAMPS  with  second- 
order  horizontal  smoothing,  (b)  ARPS  with  fourth-order  horizontal 
smoothing  coefficient  doubled  relative  to  the  control,  and  (c) 
COAMPS  with  fourth-order  accurate  horizontal  advection.  The  con¬ 
tour  interval  is  8  K. 


Fig.  7.  The  3-h  potential  temperature  for  COAMPS  using  (a)  Ax 
=  2  km  and  125  vertical  levels  and  (b)  Ax  =  2  km  and  45  vertical 
levels.  The  contour  interval  is  8  K. 


accuracy  of  the  horizontal-spatial  advection  from  sec¬ 
ond-order  (Fig.  3b)  to  fourth-order  (Fig.  6c)  increases 
the  small-scale  wave  activity  associated  with  wave 
breaking.  However,  the  fundamental  amplitude  and  ver¬ 
tical  wavelength  of  the  wave  breaking  regions  are  sim¬ 
ilar  and  appear  to  be  independent  of  hnite-difference 
technique.  Additional  experiments  indicate  a  sensitivity 
of  the  hydraulic  jump  location  and  structure  to  the  lateral 
boundary  condition  formulation  (not  shown). 

Several  experiments  were  conducted  with  degraded 
horizontal  and  vertical  resolution.  A  COAMPS  simu¬ 
lation  with  Ax  =  2  km,  shown  in  Fig.  7a,  indicates 
relative  insensitivity  to  a  doubling  of  the  horizontal  grid 
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increment.  Relative  to  the  control  (Fig.  3b),  notable  deg¬ 
radation  of  the  gravity  wave  related  structures  exists 
near  the  hydraulic  jump  and  in  the  upper-level  gravity 
wave  breaking  due  to  the  coarser  resolution.  The  3-h 
potential  temperature  obtained  from  a  COAMPS  sim¬ 
ulation  using  a  vertically  staggered  grid  with  45  levels 
and  Ax  =  2  km  is  shown  in  Fig.  7b.  The  vertical  grid 
spacing  is  based  on  what  a  typical  high-resolution  op¬ 
erational  model  might  use.  The  resolution  in  the  lower 
stratosphere,  Az  =  700  m,  is  clearly  insufficient  to  re¬ 
solve  the  wavebreaking  structures  relative  to  the  control 
simulation  (Fig.  7a),  especially  considering  the  vertical 
wavelength  is  —3-4  km.  Clearly,  some  caution  should 
be  exercised  when  choosing  the  vertical  grid  conhgu- 
ration  for  numerical  models  used  to  predict  the  occur¬ 
rence  of  wave  breaking  and  clear-air  turbulence. 

The  sensitivity  of  orographically  generated  wave 
breaking  to  the  initial  state  was  tested  by  performing 
two  simulations.  The  original  Grand  Junction  sounding 
(Fig.  1)  contains  small  vertical  wavelength  u,  t,  and  R- 
structures  in  the  lower  stratosphere  that  may  have  an 
important  influence  on  the  upper-level  wave  breaking. 
To  examine  the  relevance  of  these  vertical  structures, 
the  first  experiment  makes  use  of  smoothed  stratospheric 
velocity  and  thermal  profiles  (dotted  lines  in  Figs.  la,b). 
The  profiles  are  modified  such  that  the  layer-mean  zonal 
momentum  and  vertical  thermal  gradient  are  preserved 
from  the  original  sounding.  The  3-h  potential  temper¬ 
ature  held  using  COAMPS  (Fig.  8a)  indicates  that  the 
vertical  wavelength  of  the  upper-level  wave  breaking 
(—3  km)  and  the  lowest  altitude  of  wave  breaking  over 
the  mountain  peak  and  lee  slope  are  similar  in  both  the 
modified  initial  condition  and  control  simulations  (Fig. 
3b).  However,  there  are  significant  differences  apparent 
in  the  local  wave-breaking  regions  in  the  stratosphere. 
For  example,  in  the  10-15  km  layer,  40  km  downstream 
of  the  mountain  peak,  the  wave  breaking  is  much  more 
vigorous  and  extends  to  a  lower  altitude  in  the  control 
than  the  simulations  that  used  smoothed  profiles.  The 
gravity  wave-breaking  response  is  spatially  more  ex¬ 
tensive  in  the  simulation  that  used  smoothed  profiles, 
particularly  in  the  upstream.  Additionally,  the  maximum 
leeside  winds  at  3  h  are  —32  m  s  '  using  the  smoothed 
sounding,  while  in  the  control  experiment  has  a  65  m 
s  '  maximum.  Overall,  the  fundamental  character  of  the 
stratospheric  breaking  is  similar,  however,  the  flnescale 
structure  of  the  wave  breaking  including  the  horizontal 
and  vertical  position  and  magnitude  of  local  R,  minima 
appear  to  be  sensitive  to  initial  condition  variations  in 
the  stratosphere.  This  result  is  in  agreement  with  Blu- 
men  (1988),  who  examined  Long’s  solution  for  a  non¬ 
linear  hydrostatic  wave  with  periodic  upstream  flow  and 
found  that  the  basic  processes  of  wave  breaking  were 
unaffected.  However,  the  level  at  which  breaking  first 
occurs  and  the  wave  displacement  were  altered  by  the 
nonuniform  upstream  flow.  Thus,  the  relatively  coarse 
vertical  resolution  commonly  used  in  both  analysis  and 
model  systems  in  the  stratosphere  may  limit  the  pre- 


Distance  =  220.0  km 


Fig.  8.  The  3-h  potential  temperature  for  COAMPS  initialized  using 
(a)  smoothed  wind  and  stratihcation  prohles  in  the  stratosphere  (dot¬ 
ted  lines  in  Figs.  la,b)  and  (b)  uniform  stratospheric  conditions  (Figs. 
lc,d).  The  contour  interval  is  8  K. 

dictability  of  wave-breaking  events  through  inadequate 
specification  of  the  upstream  initial  conditions. 

The  initial  conditions  for  the  second  sensitivity  ex¬ 
periment  make  use  of  uniform  winds  and  stratification 
in  the  stratosphere,  which  effectively  removes  the  mean 
state  critical  level  near  21  km  and  is  similar  to  the  ve¬ 
locity  profile  applied  in  the  previous  simulations  of  the 
1972  Boulder  windstorm  (e.g.,  see  Figs.  lc,d).  The  sim¬ 
ulated  potential  temperature  using  COAMPS  after  3  h 
is  shown  in  Fig.  8b.  Although  some  general  similarities 
exist  between  this  simulation  and  the  control  case  (Fig. 
3b)  with  regard  to  characteristics  of  the  downslope 
windstorm  and  low-level  hydraulic  jump,  some  striking 
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Fig.  9.  West-east  oriented  vertical  cross  section  of  potential  tem¬ 
perature  after  3  h  for  a  three-dimensional  simulation  with  (a)  free- 
slip  bottom  boundary  condition  and  (b)  surface  friction.  The  contour 
interval  is  8  K.  The  heavy  solid  lines  represent  the  1.0  Richardson 
number  contour  in  the  stratosphere. 

differences  are  apparent.  Most  notably,  the  wave  break¬ 
ing  in  the  lower  stratosphere  is  more  vigorous  than  in 
the  control  experiment,  especially  in  the  layer  above  18 
km.  The  presence  of  the  critical  layer  in  the  control 
simulation  clearly  limits  the  vertical  propagation  of 
gravity  waves  and  ultimately  acts  as  an  upper  bound  on 
wave  breaking.  Once  again,  the  local  properties  of  the 
stratospheric  wave  breaking  are  sensitive  to  the  up¬ 
stream  conditions.  Additionally,  the  windstorm  is  ~10 
m  s  ‘  weaker  using  the  smoothed  stratospheric  initial 
state  of  Figs.  lc,d  relative  to  the  model  solution  that 
uses  the  Grand  Junction  sounding. 


Fig.  10.  Vertical  cross  section  of  potential  temperature  after  3  h 
for  a  two-dimensional  simulation  with  free-slip  bottom  boundary  con¬ 
dition  and  topography  representative  of  a  transect  across  the  Conti¬ 
nental  Divide  and  Front  Range.  The  contour  interval  is  8  K.  The 
heavy  solid  lines  represent  the  1 .0  Richardson  number  contour  in  the 
stratosphere. 

Several  simulations  in  three-spatial  dimensions  were 
performed  using  COAMPS  to  better  understand  the  re¬ 
lationship  between  the  two-dimensional  simulations  of 
wave  breaking  described  above  and  wave  breaking  in 
nature.  In  three-dimensions,  the  flow  is  able  to  split 
around  or  pass  over  the  topography,  in  contrast  to  two- 
dimensional  flow  that  is  either  blocked  or  is  forced  over 
the  obstacle  (e.g.,  Peng  et  al.  1995).  Clark  and  Farley 
(1984)  compared  two-  and  three-dimensional  simula¬ 
tions  of  the  1972  Boulder  windstorm  using  a  two-di¬ 
mensional  topographic  ridge.  They  compared  the  sur¬ 
face  wave  drag  in  the  two  simulations  and  found  smaller 
values  for  the  three-dimensional  solution.  In  three  spa¬ 
tial  dimensions,  turbulent  eddies  developed  in  the  con- 
vectively  unstable  region  of  the  breaking  wave  and  re¬ 
sulted  in  a  gustiness  signature  in  the  surface  wind  speed. 
The  turbulent  energy  cascade  is  not  represented  properly 
in  two  dimensions.  It  follows  that  the  three-dimensional 
solution  should  more  accurately  capture  the  turbulent 
processes  and  presumably  the  wave  breaking.  In  the 
present  study,  the  sensitivity  experiments  in  three  spatial 
dimensions  made  use  of  a  high-resolution  topography 
dataset.  The  horizontal  grid  increment  was  2  km  with 
75  vertical  levels,  distributed  in  a  vertically  stretched 
configuration  that  ranged  from  10  m  at  the  surface  to 
400  m  in  the  lower  stratosphere.  The  model  was  applied 
without  moist  processes  and  rotation.  The  Grand  Junc¬ 
tion,  Colorado,  sounding  described  above  and  applied 
in  the  two-dimensional  simulations  was  used  to  initialize 
the  model.  A  vertical  cross  section  of  potential  tem¬ 
perature  constructed  across  the  continental  divide  for 
the  3-h  simulation  time  is  shown  in  Fig.  9a.  The  wave 
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breaking  is  considerably  reduced  in  amplitude  relative 
to  the  inviscid,  two-dimensional  results  with  topography 
representative  of  a  transect  across  the  Continental  Di¬ 
vide  and  Front  Range  (cross  section  topography  of  Fig. 
9),  shown  in  Fig.  10.  The  vertical  wavelength  and  ver¬ 
tical  position  of  the  overturning  regions  are  qualitatively 
similar.  The  regions  in  the  stratosphere  characterized  by 
R-  less  than  unity  are  diminished  in  the  three-dimen¬ 
sional  simulation  as  well.  The  three-dimensionality  acts 
to  decrease  the  overall  upper-level  wave-breaking  ac¬ 
tivity.  The  wave  breaking  is  further  diminished  in  the 
presence  of  surface  friction  in  agreement  with  the  results 
of  Richard  et  al.  (1989)  and  Olafsson  and  Bougeault 
(1997),  as  illustrated  by  the  vertical  cross  section  of  the 
three-dimensional  simulation  with  surface  friction  in¬ 
cluded  shown  in  Fig.  9b.  Surface  friction  diminishes  the 
shooting  flow  in  the  lee  and  restricts  the  hydraulic  jump 
movement  resulting  in  a  position  closer  to  the  mountain, 
which  is  in  agreement  with  the  aircraft  observations 
(Fig.  2).  The  wind  speed  is  reduced  in  a  shallow  layer 
above  the  surface  and  subsequently  results  in  damping 
the  gravity  wave-breaking  response. 

4.  Summary  and  conclusions 

A  series  of  modeling  experiments  have  been  con¬ 
ducted  as  part  of  the  pre-MAP  research  effort.  Eleven 
nonhydrostatic  mesoscale  models  were  used  in  a  two- 
dimensional  mode  to  simulate  the  11  January  1972 
Boulder  windstorm  event.  The  model  configuration  was 
standardized  as  much  as  possible  with  a  horizontal  grid 
increment  of  1  km  and  a  constant  vertical  grid  spacing 
of  200  m.  The  initial  conditions,  which  were  based  on 
the  upstream  sounding  obtained  from  Grand  Junction, 
Colorado,  contained  a  critical  level  near  21  km.  The 
complete  upstream  sounding  was  retained  in  these  ex¬ 
periments,  which  presumably  enabled  for  a  more  real¬ 
istic  simulation  of  the  stratospheric  gravity  wave  break¬ 
ing.  It  follows  that  the  simulations  using  this  sounding 
represent  an  opportunity  for  a  fresh  perspective  on  the 
1972  Boulder  windstorm  event.  Simulation  of  this  case 
is  ideal  for  testing  and  evaluation  of  mesoscale  numer¬ 
ical  model  codes  and  numerical  algorithms  because  of 
the  complex  wave-breaking  response  and  relative  in¬ 
sensitivity  to  the  upper-boundary  condition. 

The  results  from  the  simulations  are  very  encouraging 
with  regard  to  gravity  wave-breaking  predictability  in 
that  all  models  produced  upper-level  wave  breaking  in 
similar  horizontal  locations  and  vertical  layers.  The  up¬ 
per-level  wave  breaking  was  most  pronounced  in  the 
13-16  km  and  18-20  km  layer  with  an  upstream  tilt 
with  height.  The  wave  breaking  extended  horizontally 
downstream  by  —80-100  km  and  —20  km  upstream. 
All  of  the  numerical  model  simulations  contained  strong 
leeside  winds  in  the  low-levels  with  a  well-defined  hy¬ 
draulic  jump,  in  agreement  with  the  observations  of 
Lilly  and  Zipser  (1972),  Klemp  and  Lilly  (1975),  and 
Lilly  (1978).  However,  significant  differences  among 


the  simulations  were  apparent  in  the  movement  of  the 
hydraulic  jump  and  the  depth  of  the  breaking  layers. 
The  low-level  wind  speed  and  vertical  velocity  maxima 
in  the  lee  varied  by  a  factor  of  2  among  the  simulations 
at  specific  integration  times.  However,  in  some  circum¬ 
stances  a  comparison  of  global  maxima  and  minima 
among  models  may  be  difficult  to  interpret  in  the  pres¬ 
ence  of  high-frequency  transients. 

Sensitivity  experiments  indicate  that  the  structure  of 
the  wave  breaking  was  impacted  by  the  numerical  dis¬ 
sipation,  numerical  representation  of  the  horizontal  ad- 
vection,  and  the  lateral  boundary  conditions.  The  upper- 
level  wave  breaking  was  found  to  be  particularly  sen¬ 
sitive  to  the  vertical  resolution,  implying  that  the  vertical 
grid  spacing  typically  used  in  operational  models  is  like¬ 
ly  insufficient  to  resolve  most  upper-level  gravity  wave 
propagation  and  breaking  processes.  Substantial  differ¬ 
ences  exist  in  the  simulation  of  the  wave  breaking  be¬ 
tween  one  of  the  models,  which  did  not  parameterize 
the  subgridscale  eddy  viscosity,  and  the  others  models, 
which  made  use  of  turbulence  parameterization 
schemes.  The  simulations  of  upper-level  wave  breaking 
were  sensitive  to  the  upper-level  upstream  conditions. 
Variations  in  the  vertical  shear  and  stability  above  10 
km  contributed  to  significant  changes  in  the  structures 
associated  with  wave  breaking.  Additionally,  three-di¬ 
mensionality  and  surface  friction  act  to  decrease  the 
amplitude  of  the  upper-level  breaking. 

The  results  from  multimodel  simulations  suggest  that 
the  turbulent  breakdown  of  orographically  generated 
gravity  wave  episodes  may  be  predictable  in  some  cir¬ 
cumstances.  Wave  breaking  is  generally  thought  to  be 
a  highly  nonlinear  phenomenon  (e.g.,  Peltier  and  Clark 
1979)  and  a  challenge  for  numerical  models  to  accu¬ 
rately  simulate.  However,  the  coherent  transients  and 
longer-temporal  scale  structures  within  the  gravity  wave 
breakdown  regions  are  surprisingly  predictable  as 
evinced  by  the  similarities  among  the  numerical  simu¬ 
lations.  The  results  from  this  test  case  are  based  on  a 
group  of  two-dimensional  highly  idealized  models.  The 
particular  case  chosen  is  well  suited  to  two-dimensional 
idealization.  Three-dimensional  wave  breaking  in  the 
presence  of  temporally  varying  complex  forcing  is  a 
significantly  more  complicated  environment  to  evaluate 
wave-breaking  predictability.  Three-dimensional  test 
cases  are  currently  being  performed  that  incorporate 
complex  topographical  forcing  such  as  the  Alps,  as  well 
as  real-data  initial  states,  moist  dynamics,  and  rotation. 
With  these  numerical  experiments  performed  in  con¬ 
junction  with  the  MAP  research  effort,  confidence  may 
be  established  in  the  capability  of  sophisticated  non¬ 
hydrostatic  models  to  predict  episodes  of  the  turbulent 
breakdown  of  topographically  generated  gravity  waves. 
These  models  will  serve  to  guide  research  missions  to 
observe  wave  breaking  during  the  field  phase  of  MAP 
and  will  serve  as  dynamical  tools  to  increase  our  ru¬ 
dimentary  understanding  of  gravity  wave  dynamics  over 
three-dimensional  complex  topography. 
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